Self-consistent models of triaxial galaxies in MOND gravity 

You-gang Wang^'^, Xufen Wu^, HongSheng Zhao^'^ 
hz4(§st-andrews .ac.uk 

ABSTRACT 

The Bekenstein-Milgrom gravity theory with a modified Poisson equation is 
tested here for the existence of triaxial equihbrium solutions. Using the non- 
negative least square method, we show that self-consistent triaxial galaxies exist 
for baryonic models with a mild density cusp p ~ ^. Self-consistency is achieved 
for a wide range of central concentrations, S ~ 10 — IOOOMqPc"^, representing 
low-to-high surface brightness galaxies. Our results demonstrate for the first time 
that the orbit superposition technique is fruitful for constructing galaxy models 
beyond Newtonian gravity, and triaxial cuspy galaxies might exist without the 
help of Cold dark Matter. 

Subject headings: galaxies: kinematics and dynamics- methods: numerical - dark 
matter 



1. Introduction 

Constructing models of galaxies in triaxial equilibrium is a classical challenge in dy- 
namics, either in Newtonian or Modified gravity. However, due to numerical construc- 
tions, few such models have been developed in Newtonian since the pioneering work of 
Schwarzchild(1979) on triaxial elliptical galaxies, and Zhao (1996) on the fast rotating tri- 
axial bar of the Milky Way. Merritt & Fridman (1996) (hereafter MF96) considered a model 
with a central density cusp and found that it is possible to construct a model of triaxial 
equilibrium galaxies self-consistently in Newtonian dynamics. Models with a steeper density 
cusp are not in equilibrium. This has generated interests in understanding cusp-triaxiality 
relation, the triaxiality-velocity anisotropy relation and the cusp-black hole relation. Re- 
cently, Capuzzo-Dolcetta et al. (2007) modeled cuspy triaxial galaxies in A Cold Dark 

^National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100012 
2 SUPA, School of Physics and Astronomy, University of St Andrews, KY16 9SS, UK 
■^Department of Astronomy Peking University, Beijing 100871, China 



- 2 - 



Matter (CDM) haloes. They took the model of MF96 as the luminous density distribution, 
and found that a model of cupsy triaxial galaxies with CDM halos is also self-consistent. 

A tough problem for ACDM (A plus Cold Dark Matter) is that there is no physics 
to justify A the tiny cosmological constant or dark energy (see White 2007, Sarkar 2007). 
Although future particle physics experiments might well prove the existence of new species 
of particles and vacuum energy, their experimental-determined abundance could easily be 
factors of a few different from the 1:3 ratio as precisely required by fitting observations 
of Microwave background (Zhao 2006). This motivates exploring alternative theories like 
the CO- variant version of Modified Newtonian Dynamics (MOND) in the mean-time. E.g., 
instead of more traditional interpretation of modified gravity (Bekenstein 2004), MOND was 
given the interpretation as a co- variant Dark Energy (DE) the VA model of Zhao (2007). 
This model uses General Relativity plus a non-uniform Dark Energy (DE) fluid (described 
by a four- vector flow) to give both the effects of DE and DM. Perturbations of such DE fluid 
bends space-time as an effective DM without actually invoking DM. Unlike ACDM halos, 
which are positive and nearly round, the effective DM can sometimes be extremely flattened 
or negative in some regions of galaxies or clusters (Wu et al. 2007, 2008), although the data 
are not good enough to tell such negative regions yet (Nipoti et al. 2007a). 

Overall MOND and ACDM halos are comparably successful in explaining the flat rota- 
tion curves of high surface brightness discs and satellites orbits at large radii (e.g., Angus et 
al. 2008). However, both theories are faced with severe challenges: ACDM hinges on uncer- 
tain feedback on galactic scales, and even the maximum feedback falls short of explaining 
the velocity curves of low-surface brightness systems (Gnedin & Zhao 2002); Bekenstein's 
(2004) co-variant MOND seems to rely on non-relativistic neutrinos on large scale to match 
the weak lensing data (Angus et al. 2007) and the Cosmic Background Radiation (Skordis 
et al. 2006). 

The tight correlation observed between the mass profiles of baryonic matter and dark 
matter at all radii in spiral galaxies (e.g., McGaugh et al. 2007; Famaey et al. 2007) is 
still the best case for the MOND paradigm of Milgrom (1983), which postulates that for 
accelerations below = 1.2 x lO^^'^ms"^ the effective gravitational attraction approaches 
idNO-oY^'^, where is the usual Newtonian gravitational field. Indeed, without resorting to 
galactic dark matter, this simple prescription is amazingly successful at reproducing galactic 
rotation curves over five decades in mass ranging from tiny dwarfs (e.g.. Gentile et al. 2007) to 
early-type disk galaxies (e.g., Sanders & Noordermeer 2007) to massive ellipticals (Milgrom 
& Sanders 2003). Some outliers exist in models of gravitational lensing (Zhao et al. 2006, 
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Chen & Zhao 2006, Shan et al. 2008). Lensing and structure formation are well-defined 
calculations in co- variant MOND (Halle & Zhao 2007). Galaxy formation simulations in 
classical MOND have been carried out too, and it is possible to form bars (Tiret & Combes 
2007), and elliptical galaxies by mergers (Nipoti et al. 2007b), and spherical bulges by 
instability (Zhao, Xu & Dobbs 2008). 

Very little is known about the triaxial stationary equilibrium in non-Newtonian dy- 
namics, although one speculates that such equilibrium might exist. In this paper, we test 
the existence of triaxial galaxies in the Bekenstein-Milgrom MOND theory, a well-defined 
classical theory with relativisitic and cosmological extentions. We apply the same density 
model given by MF96 with a fixed mass, and axis ratios. Our aim is to examine whether 
this model is self-consistent in MOND. 

The MOND theory is scale-dependent. When the gravitational acceleration is below oq, 
the scaling deviates from the Newtonian law. Here we consider models with different 
ratio of oq to the simulation unit G x (Total Mass) / (Scale Radius)^. We start with a 
nearly Newtonian system (high surface brightness) and finish with a deep-MOND system 
(low surface brightness). 

The rest of the paper is organized as follows. In §2, we present the density distribution of 
the galaxy and the corresponding potential and forces. The orbits are calculated in §3. The 
construction of self-consistent models is discussed in §4, and in §5, we give our conclusions 
and disscussions. 



2. The density, potential distribution and forces in the systems of triaxial 

galaxies 

The baryon component adopted by Capuzzo-Dolcetta (2007) has the following ellipsoidal 
mass distribution: 

Ml 
^'^^ ~ 2'Kh,hyK R{1 + Rf 

where -R^ = + -l- ^{hz ^ hy < h^) is the length of the major axis, and M is the total 
baryon mass. For the ratio of the three major axes, we adopt : hy : = 1 : 0.86 : 0.7 
throughout this paper. This density profile was studied in MF96, but in their work, it 
represented the density profile of the dark matter. The density profile is plotted in the top 
panel of Fig. [TJ Density is cuspy and scales as at the center and decreases quickly as 
i?""^ at large radii. 



Hereafter, we adopt units in which the total mass M, the x-axis scale length hx, and 
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the gravitational constant G are unity. Thus one unit of time corresponds to 

. / M \-V2/ n3/2 

= 1.49xl0V(^i ) {—^] 2 



The MOND nonrelativistic field equation is different from Poisson's equation in Newton's 
theory (Bekenstein & Milgrom 1984), and it has the following form 

AnGp = V-(;uV$) 

= d,{fid,<l>)+dy{fldy<l>)+d,{fid,<^), (3) 

where $ is the gravitational potential generated by the density distribution p. We adopt the 
function p as in Zhao & Famaey (2006), which is given by 

M = (4) 
9 + ao 

Clearly p = 1 for g ^ oo and p = g/a^ for (7 ^ 0, as expected for any MOND theory. The 
gravity strength g is given by 

9 = |V$| 



{d^^f + {dy^Y + {d^^f. (5) 

From equation ([5]) it can be seen that the equation ([3]) does not reduce to Poisson's 
equation if oq is not small enough. We solve equation ([3]) numerically, using the Bologna 
group's MOND Poisson Solver (Ciotti et al. 2006), which is based on a spherical grid. This 
code yields results which are consistent with that of the Cartesian grid code of the Paris group 
(Tiret & Combes 2007). We apply the Poisson Solver with a resolution of 256 x 64 x 128 
grids points, and choose the radial grid points as = 2.0tan[(i + 0.5)2||fi]kpc. We obtain 
the components of gravity in x, y, z directions, ie. values of d^-, dy, and dz from the Poisson 
Solver. 

In the lower panel of Figure 1, we show the potential radial profile along the major axis 
for different values of oq in simulation units. In the middle panel we show the radial gravity 
\gr\/cLo as a function of R along the long axis. In the simulations is dimensionless and is 
equal to 0.0833, 0.833 and 8.33, which correspond to simulation units of m = 1.0 x lO^^M©, 
1.0 X IO^Mq, and 1.0 x lO^M© respectively, and G = 1. In all simulations, we start with a 
total mass of M = 10m (See Tabled]). We find that there is no obvious difference between 
the mild and weak gravity cases, except in the central region. For the strong gravity case 
(ao = 0.083), the gravitational potential falls quickly at the outer regions of the galaxy. 
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The left panel of Figure [2] shows the isopotentials (solid lines) in the gas with axial 
ratios : hy : = 1 : 0.86 : 0.7 at five different radii, R =2, 10, 20, 30 and 50 kpc. There 
are no obvious differences with different values of ao(see Fig. Further more, comparing 
with the density contours (the dotted lines on the Fig. [2^), the isopotential contours are 
more spherical than the isodensity (Fig. [2]d), which is expected (see also the examples in 
Lee & Suto (2003) and Wang & Fan (2004)). The three-dimensional potential distribution 
is approximately ellipsoidal, with eccentricity increasing from the center to the outer part 
of the cluster. For example, the axial ratio of potential distribution is : h^. = 0.90 : 1.0 
at 2kpc, while : = 0.99 : 1.0 at 100 kpc. The dashed lines are isodensities of the 
effective DM halo in MOND. On the Fig. [2]3, it is interesting that the axis ratio of effective 
DM densities are smaller than that of baryon densities at radii smaller than 1 kpc, where 
baryons dominate the dynamics. To produce the same dynamical behavior(e.g. the rotation 
curves of galaxies), the CDM halos always adopt a constant axis ratio of dark matter density 
in the inner and outer part of galaxies, otherwise the CDM models would be much too 
complex. 

For the model of the elliptical galaxy, one important parameter is the effective radii. 
Hernquist (1990) shew that the effective of the elliptical galaxy is ~ l.SlbSh^ for a spher- 
ical system. In our models, the effective radii is also ~ l.Sh^ because this parameter is 
independent of gravity and is set by the baryon density profile. 

MOND is equivalent to an effective DM theory. The effective density can be obtained 

from 

47rGp^^/, = V-[(l-/x)V«l']. (6) 

This density can be quite different from both the baryonic density and also the nearly 
ellipsoidal density of the CDM halo. A example is shown in Fig. [2l 

3. Integration of orbits 

In Newton's theory, energy, and the three components of the angular momentum, are 
integrals of motion in a spherical system. However, only energy is kept as a constant in the 
triaxial potential. MOND does not change the number of constants of motion. Bekenstein- 
Milgrom's theory conserves energy for a stationary potential. In order to obtain a specified 
accuracy, we use the 7/8 order Runge-Kutta algorithm (Fehlberg 1968) to carry out the 
integrations. The integration time for a full orbit is 100 dynamical times. One dynamical 
time is defined as the period of the 1:1 resonant orbit in the x-y plane. Table [2] lists the 
periods of 1:1 resonant orbits in the x-y plane as a function of energy. 
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In order to produce a library of all possible orbits, we followed Schwarzschild (1993) 
and MF96 in selecting our initial conditions. Specifically, the definition of start-space is the 
same as in MF96. The orbital catalogues include both initially stationary orbits and x-z 
plane launched orbits with nonzero velocity Vy, where Vy = \/2{E — $) 7^ is the velocity 
along the y-axis. For each model we have 20 energy surfaces. The surfaces are the critical 
surface of the shells which divide the total system into 21 parts of equal mass. One energy 
consists 342 initial conditions, where 192 orbits are from the stationary start-space and the 
rest are from the x-z start-space. Therefore, for each model we have 6840 orbits. 

Figure [3] shows the orbits in three separate planes for MOND models. The energies of 
orbits are from shell 10. The orbits in a triaxial MOND potential are similar to the orbits 
in a triaixal Newtonian potential, which yield the box, tube and stochastic orbits. 

4. Construction of self-consistent models 

Since Schwarzschild (1979) proposed a numerical method to construct self-consistent 
models of galaxies, orbit superposition techniques have become an important tool in dynam- 
ical modelling (e.g., Rix 1997; van der Marel et al. 1998; Kuijken 2004; Binney 2005; van de 
Yen 2006; Capuzzo-Dolcetta et al. 2006). The procedure of constructing such models can 
be described as follows: 

(1) A three-dimensional mass distribution is employed and the corresponding gravita- 
tional potential is calculated by solving Poisson's equation. In this paper, the case is revised 
as we solve the modified Poisson's equation. The simulations are also not scale-free, which 
was possible for some Newtonian galaxy models. 

(2) A full orbital library is established. In other words, it is necessary to generate a large 
set of initial conditions for the orbits. 

(3) The whole system is divided into many cells with the same mass. The fraction of 
time spent by each orbit in each cell Oij is recorded into a large matrix. 

(4) The non-negative occupation numbers Wj are programmed in a linear equation group 
with the above matrix and the desired mass in each cell. The weights are then solved by 
non-negative least squares inversion. Note that the relation between orbital weights and the 
density distribution is linear independent of gravity theories. 

The linear equation group can be written as 
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No 

J2WjOij = M„ 1 = 1, ,iVe, (7) 

where Mi is the mass of each grid cell, No is the number of orbits, here we have A'";, = 6840. 
Nc is the number of cells, and we adopt A^c=912. Wj > is the weight number of the 
occupation. We follow MF96 to divide the first octant into 960 cells, which can be described 
by the following steps. The galaxy was divided by 20 ellipsoidal shells into 21 sections of equal 
mass. For each shell, only the first octant was considered. The first octant was then divided 
into three parts by the planes z = hzx/hx, y = hyx/h^, and z = hzv/hy. Each part, for 
example the first part, which contained the x-axis, was divided by h^y/hyX = 1/5,2/5,2/3 
and h^z/h^x = 1/5,2/5,2/3 into 16 cells. The total number of cells is 20 x 3 x 16 = 960. 
We consider the inner 912 cells in solving the linear equation groups. 

Generally, equation ([7]) can be solved by linear programming (e.g., Schwarzschild 1979, 
1982, 1993) using Lucy's method (Lucy 1974; Statler 1987), or maximum entropy methods 
(e.g., Richstone & Tremaine 1988; Statler 1991; Gebhardt et al. 2003), or with the least 
squares solver (Lawson & Hanson 1974) (MF96; van de Yen et al. 2006; Capuzzo-Dolcetta 
et al. 2006). The four methods have their own advantages and disadvantages, which we do 
not discuss here. We adopt the least squares method to solve the linear programming. Then 
the minimum square deviation can be written as 

.AT, No 2 

X' = J^J2{C^-J2^,0,) (8) 

" j=l i=l 

For the symmetry of the triaxial model, only one octant Oij is considered in solving the 
linear programming. 



5. Conclusion and discussion 

In practice, a parameter 6 is adopted to describe the departure form self-consistency, 
which is defined as in MF96: 



where M is the average mass in each cell. Figure H] shows the departure from self-consistency 
as a function of the number of orbits. Note that there is a remarkably strong dependence on 
the number of orbits. The departure parameter S drops quickly when the number of orbits 
exceeds 4000. 6 is 10~^^ when the 6840 orbits are adopted in the optimization routine. 
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which shows that all three triaxial models are self-consistent in MOND. Amazingly, we also 
find that the departure parameter 6 is insensitive to ag. In other words, if one of the triaxial 
galaxies is self-consistent in MOND, then the same model with a different value of will 
also be self-consistent. 

Figure [5] shows the contribution of mass for various orbital families in the self-consistent 
model. Chaotic orbits contribute a larger fraction of mass with the increase of energy. This 
shows a certain degree of consistency with that of Capuzzo-Dolcetta et al. (2007). However, 
discrepancies between our results and that of Capuzzo-Dolcetta et al. (2007) are apparent. 
There is a clear variational trend of the cumulative energy distribution of the various orbital 
families in Newton's case, however, no apparent variational trend for different orbits has 
been found in MOND. 

In Figureini we show the radial profile of the rescaled radial velocity dispersion cr,./Kir,oo(left 
panel) and the radial profile of the anisotropy parameter P = 1 — 0.5(< V"^ — V^ >)/ < Vj^ > 
(right panel), where V^ir,oo is the circular velocity at infinity, and Vr is radial velocity. The 
parameter V"^ is defined as = + Vy + v"^, where v^, Vy and are the three components 
of velocity. It can be seen that the profiles of the radial velocity dispersion are nearly fiat 
in MOND, except in the central region. Recently, Capuzzo-Dolcetta et al. (2007) studied 
the same model as that used here, but in the Newtonian gravity. It can be seen that the 
profile of the velocity dispersion in our model shows a certain degree of consistency with 
that of Capuzzo-Dolcetta et al. (2007), but the discrepancies between our results and that 
of Capuzzo-Dolcetta et al. (2007) are apparent. There is central peak in our velocity dis- 
persion profile. The global average of the velocity is 170.49, 78.84 and 42.70 km/s for 
ao=0.083, 0.833, and 8.333, respectively (see the last column of Table [1]). 

The anisotropy parameter (3 > suggests that the stars have high radial velocity and 
the ratio of < > / < — > reaches a maximum value at 2 — 3 kpc. Our anisotropy 
profiles are similar to that of Nipoti et al. (2007b) form dissapationless collapse in MOND 
(their figur 2) and in Newtonian gravity van Albada (1982) and Wang et al. (2008a). A 
gently rising beta is widely used (e.g., Milgrom & Sanders 2003). However, the peak in the 
anisotropy profile in our model is not seen in any simulations and observations to our best 
knowledge. This is might be indication of problem of MOND. Our models have a sudden 
change of the amount of box orbit near shells 10-12 (radius 2-3kpc), which could be the 
reason. 

In order to check the stability of our models, we used the Antonov's third law (Binney 
& Tremaine 1987). If the stellar system is stable, then its density p and potential $ satisfy 
everywhere the inequality d^pb/d$^ < 0. Figure [7] shows that d^pb{0,y,0)/d\^{0,y,0)\^ is 
positive everywhere along the intermediate axis for our three models, since the values of 
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potential are negative, our models are consistent with the law, which means the models 
appear globally stable. We will also be checking the stability of the models in collaboration 
with groups which have MOND N-body code. Effects due to pattern rotation and external 
field in MOND remain to be studied (Wang et al. 2008b). 

We have used the Schwarzschild approach to examine whether the triaxial galaxies 
are self-consistent in MOND. Using the Bologna Poisson solver to determine the potential 
and the accelerations at some discrete points, we used three-dimensional interpolation then 
obtain the potential and accelerations at arbitrary point. The orbits conserve energy to a 
very high level of accuracy, which shows that our method is feasible in calculating orbits 
using a grid-based MONDian potential. The departure parameter 6 shows that the triaxial 
galaxy models adopted in this paper in MOND are self-consistent. Our results show that it 
is theoretically allowed to have triaxial galaxies in rigorous equilibrium in a universe with 
either non- Newtonian gravity or certain models of Dark Energy (Zhao 2007). 
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Table 1: Main physical parameters for our three ellipsoidal MOND models 



Model 
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Fig. 1. — Upper: the radial distribution of the dimensionless density p in units of 
M/ {2'n'hxhyhz), which is identical for the three models. Middle: the MOND gravity strength 
|(?r|/ao as a function of R along the long axis. Lower: the potential major axis profile, for 
different values of Oq, obtained from the Poisson solver. 
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Fig. 2. — Left panels: Isopotentials (solid lines), isodensities of effective dark matter (dashed) 
and isodensities of baryons (dotted) on the x-z plane, in the model of = 0.833. The con- 
tours of the other two models are almost same. The isopotentials are essentially a snapshot 
of the metric perturbations, and the isodensities of effective DM halos are images of the 
perturbation of certain Dark Energy fluid (Zhao 2007). Right panel shows the axis ratio 
c/a (the ratio between the short and long axis) of the potentials (upper series of curves), 
effective DM densities (lower series of curves) and baryon densities (the long-dashed hne) in 
a log-log space. For each series, the sohd, dotted and dashed lines are for different models 
of oo as in the legend. The oscillations at large radii are due to the numerical resolution. 
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Fig. 3. — Orbits in a nonrotating triaxial potential in MOND. The left panels are for a 
stationary start space, while the right panels are for a x-z start space (1:1 resonant orbit). 
The energy is from the 10th cell. 
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Fig. 4. — Departure form self-consistency 5 as a function of the number of orbits for different 
values of oq. The triangles, asterisks, and diamonds represent oq = 0.083, 0.833, and 8.333, 
respectively. 
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Table 2: The energy and dynamical time in every shell 
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Fig. 5. — Orbital families in the self-consistent models. The symbols b, x, z, and c denote 
the mass contributed by box, x-tube, z-tube, and chaotic orbits, respectively. 
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Fig. 6. — The radial profiles of the rescaled radial velocity dispersion o"r./K;ir,oo (left panel) 
and of the anisotropy parameter (3 = 1 — 0.5 < V"^ — > / < > (right panel). The 
solid, dotted and dashed lines represent Oq = 0.083, 0.833, 8.333, respectively. 
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